#!usr/bin/perl -w

use strict;

die "fa.list format\nxx1\txx1.fa\tcorr-name1\nxx2\txx2.fa\tcorr-name2\n..\nusage:perl $0 <fa.list> <all_orthomcl.out> \n" unless @ARGV == 2;

open A,$ARGV[0];    open B,$ARGV[1];	open FAM,">gene.family.xls";	print FAM "Orth_id\tfam_id";
my $num=0;my %h; my %gene; my %gene_id;  my %res;my %core;	my @aa; my @sam;   my %num;my %map;	my %use;	my %handle1; my %handle2;    my %handle3; 
my %gene_id2;	my %map_name;
my $shelln="special_gene";  mkdir $shelln unless (-d $shelln);	my $homolog=0;
my %hash2;
while(<A>){
	chomp;my @a=split;  $num++;  	open IN,$a[1];	my $temp=0;	$map_name{$a[0]}=$a[0];
	$/=">";    <IN>;	
	if(@a == 3){
		push @sam,$a[0];    print FAM "\t$a[2]";				 $map_name{$a[0]}=$a[2];
		my $handle2;    open $handle2,">$a[2].imp";              $handle2{$a[0]}=$handle2;
		my $handle3;    open $handle3,">$shelln/$a[2].special.fa";       $handle3{$a[0]}=$handle3;
	}
	else{
		push @sam,$a[0];    print FAM "\t$a[0]";				 $map_name{$a[0]}=$a[0];
		my $handle2;    open $handle2,">$a[0].imp";              $handle2{$a[0]}=$handle2;
		my $handle3;    open $handle3,">$shelln/$a[0].special.fa";       $handle3{$a[0]}=$handle3;
	}
	$core{$map_name{$a[0]}}[1]=0;	$core{$map_name{$a[0]}}[0]=0;
	while(my $inf=<IN>){
		$temp++;chomp $inf;my @b=split /\n+/,$inf;my @c=split /\s+/,$b[0];
		if(@a == 3){
			$gene_id{$c[0]}=$a[2];	
		}
		else{
			$gene_id{$c[0]}=$a[0];
		}

		$h{$c[0]}=$inf;	$gene_id2{$a[0]}{$c[0]}=1;
		$core{$map_name{$a[0]}}[1]++;
	}
	close IN;
	#$num{$a[0]}=$temp;
	$/="\n";
}
close A;	print FAM "\n";
my $rec=0;  my $shell1="single.copy.gene";  mkdir $shell1 unless (-d $shell1);	my $order=-1;
open OUTT,">single.copy.genefamily.xls";    open KK,">all.orth.gene.pro.fa"; open MUS,">muscle.list";
while(<B>){
	chomp;my @a=split /\t+/;	my @b=split /[(, ]/,$a[0];	my @c=split /\s+/,$a[-1];	$order++;
	if($b[3] == $num){
		$homolog++;
	}
	if($b[1] == $num && $b[3] == $num){
		$rec++;	print OUTT "$b[0]";
		open OUT,">$shell1/$b[0].fa";
		foreach my $k(1..$#c){
			my @d=split /[\(\)]/,$c[$k];
			print OUT ">$h{$d[0]}";			print OUTT "\t$map_name{$d[1]}|$d[0]";
		}
		close OUT;
		`muscle -in $shell1/$b[0].fa -out $shell1/$b[0].fa.out -clw -quiet`;
		open  PART,"$shell1/$b[0].fa.out";		open  OUT_P,">$shell1/$b[0].out.fa";		print MUS "$b[0]\t$shell1/$b[0].out.fa\n";
		my %part=();<PART>;	<PART>;
		while(my $inf=<PART>){
			chomp $inf;			next if($inf=~/^ /);	next if($inf=~/^\s/); 	next if($inf=~/^\t/);
			my $len=length $inf;
			next if($len < 2);
			my @dd=split /\s+/,$inf;
			#next if(@dd != 2);
			$part{$dd[0]}.=$dd[-1];
		}
		close PART;
		foreach my $k1(sort keys %part){
			$res{$gene_id{$k1}}.=$part{$k1};			print OUT_P ">$k1\n";            print OUT_P "$part{$k1}\n";
		}
		close OUT_P;
		print OUTT "\n";
	}
	if($b[3] == 1){
		foreach my $k(1..$#c){
			my @d=split /[\(\)]/,$c[$k];			my $out3=$handle3{$d[1]};			print $out3 ">$h{$d[0]}";	$core{$map_name{$d[1]}}[0]+=$b[1];
		}
	}
	my %fam;	my $subrec=0;	
	foreach my $k1(0..$#c){
		my @d=split /[\(\)]/,$c[$k1];next if(@d < 1);
		$map{$order}{$d[1]}++;		$gene{$d[1]}{$order}{$d[0]}++;		$fam{$d[1]}++;	
	}
	for (my $i=0;$i< $#sam;$i++){
		for (my $j=$i+1;$j<=$#sam;$j++){
			if(exists $fam{$sam[$i]} && exists $fam{$sam[$j]}){
				$hash2{$i}{$j}++;
			}
		}
	}
	print FAM "$b[0]\t$order";	
	foreach my $k1(@sam){
		if(exists $fam{$k1}){
			print FAM "\t$fam{$k1}";	$subrec++;
		}
		else{
			print FAM "\t0";
		}
	}
	print FAM "\t$_\n";
}
close B;close OUTT; close MUS;
`rm single.copy.gene/*.fa.out`;
foreach my $k(sort keys %res){
	my $len=length $res{$k};
	print KK ">$k\n";
	for(my $i=0;$i <= $len;$i+=60){
		my $subseq=substr($res{$k},$i,60);
		print KK "$subseq\n";
	}
}
close KK;
#`muscle -maketree -in all.orthology.gene.pro.merge.fa -out all.orthology.gene.pro.merge.fa.tre -cluster neighborjoining`;
if(@sam <= 5){
	my $sam=join ",",@sam;      my $sam2=join ".imp,",@sam;
	print STDERR "perl /mnt/sdb/zengl/bin/rna-pipe/venn/Venn-RNAseq.pl -f $sam2.imp -l $sam\n";
}
my %record;
# specail
foreach my $k1(sort {$a<=>$b} keys %map){
	my @kk=sort keys %{$map{$k1}}; my $name;
	foreach my $k2(@kk){
		$use{$k2}+=$map{$k1}{$k2};
		my $out2=$handle2{$k2};
		if(@kk >= 2){
			$name="tmp_".$k1;
		}
		if(@kk == 1){
			$record{$k2}++;			$name=$k2.".sam_".$record{$k2};
		}
		print $out2 "$name\n";
		foreach my $k3(sort keys %{$gene{$k2}{$k1}}){
			delete $gene_id2{$k2}{$k3};
		}
	}
}
#print special gene info
foreach my $k1(sort keys %gene_id2){
	foreach my $k2(sort keys %{$gene_id2{$k1}}){
		$record{$k1}++;
		my $out2=$handle2{$k1};			my $out3=$handle3{$k1};
		print $out2 "$k1.sam_$record{$k1}\n";		print $out3 ">$h{$k2}";
		$core{$map_name{$k1}}[0]++;
	}
}

open CORE,">gene.family.info.xls";
print CORE "Sam\tsingle.copy.gene\tspecial_gene\tall\n";
foreach my $k(sort keys %core){
	print CORE "$k\t$rec\t$core{$k}[0]\t$core{$k}[1]\n";
}
close CORE;	close FAM;
open MAT,">genefamily.strain.mat.xls";
foreach my $k1(0..$#sam){
	print MAT "$sam[$k1]";
	foreach my $k2(0..$#sam){
		if (exists $hash2{$k1}{$k2}){
			print MAT "\t$hash2{$k1}{$k2}";
		}else {
			print MAT "\t-";
		}
	}
	print MAT "\n";
}
print MAT "All.strains";
foreach my $k1(1..$#sam){
	print MAT "\t-";
}
print MAT "\t$homolog\n";

